function [f,g] = dgamma(c,x)
% DGAMMA Derivative of Gamma PDF with respect to alpha
%	[F, G] = DGAMMA(ALPHA, X)
%	Return partial derivative and second order partial derivative
%	This is an utility function for GAMMLE

     t1 = [-10.42375494, ...
    -5.28903990, ...
    -3.50252422, ...
    -2.56138454, ...
    -1.96351003, ...
    -1.54061921, ...
    -1.22002355, ...
    -0.96500857, ...
    -0.75492695, ...
    -0.57721566, ...
    -0.42375494, ...
    -0.28903990, ...
    -0.16919089, ...
    -0.06138454, ...
    0.03648997, ...
    0.12604745, ...
    0.20854787, ...
    0.28499143, ...
    0.35618416, ...
    0.42278434, ...
    0.48533597, ...
    0.54429344, ...
    0.60003988, ...
    0.65290117, ...
    0.70315664, ...
    0.75104745, ...
    0.79678317, ...
    0.84054699, ...
    0.88249995, ...
    0.92278434, ...
    0.96152644, ...
    0.99883889, ...
    1.03482249, ...
    1.06956784, ...
    1.10315664, ...
    1.13566284, ...
    1.16715354, ...
    1.19768985, ...
    1.22732754, ...
    1.25611767, ...
    1.28410709, ...
    1.31133889, ...
    1.33785279, ...
    1.36368548, ...
    1.38887093, ...
    1.41344062, ...
    1.43742381, ...
    1.46084774, ...
    1.48373779, ...
    1.50611767, ...
    1.52800953, ...
    1.54943413, ...
    1.57041093, ...
    1.59095821, ...
    1.61109315, ...
    1.63083192, ...
    1.65018977, ...
    1.66918107, ...
    1.68781943, ...
    1.70611767, ...
    1.72408796, ...
    1.74174182, ...
    1.75909018, ...
    1.77614340, ...
    1.79291133, ...
    1.80940335, ...
    1.82562836, ...
    1.84159487, ...
    1.85731095, ...
    1.87278434, ...
    1.88802239, ...
    1.90303214, ...
    1.91782034, ...
    1.93239340, ...
    1.94675748, ...
    1.96091850, ...
    1.97488209, ...
    1.98865369, ...
    2.00223849, ...
    2.01564148, ...
    2.02886746, ...
    2.04192103, ...
    2.05480664, ...
    2.06752853, ...
    2.08009082, ...
    2.09249745, ...
    2.10475222, ...
    2.11685882, ...
    2.12882077, ...
    2.14064148, ...
    2.15232425, ...
    2.16387225, ...
    2.17528856, ...
    2.18657615, ...
    2.19773788, ...
    2.20877652, ...
    2.21969475, ...
    2.23049518, ...
    2.24118032, ...
    2.25175259];

    t2 = [101.43329915, ...
    26.26737721, ...
    12.24536455, ...
    7.27535659, ...
    4.93480220, ...
    3.63620967, ...
    2.83404916, ...
    2.29947414, ...
    1.92253996, ...
    1.64493407, ...
    1.43329915, ...
    1.26737721, ...
    1.13425343, ...
    1.02535659, ...
    0.93480220, ...
    0.85843189, ...
    0.79323283, ...
    0.73697414, ...
    0.68797206, ...
    0.64493407, ...
    0.60685287, ...
    0.57293276, ...
    0.54253746, ...
    0.51515251, ...
    0.49035776, ...
    0.46780689, ...
    0.44721207, ...
    0.42833216, ...
    0.41096375, ...
    0.39493407, ...
    0.38009550, ...
    0.36632119, ...
    0.35350154, ...
    0.34154140, ...
    0.33035776, ...
    0.31987790, ...
    0.31003786, ...
    0.30078114, ...
    0.29205768, ...
    0.28382296, ...
    0.27603723, ...
    0.26866494, ...
    0.26167418, ...
    0.25503621, ...
    0.24872510, ...
    0.24271741, ...
    0.23699184, ...
    0.23152906, ...
    0.22631146, ...
    0.22132296, ...
    0.21654883, ...
    0.21197560, ...
    0.20759089, ...
    0.20338331, ...
    0.19934239, ...
    0.19545843, ...
    0.19172249, ...
    0.18812629, ...
    0.18466215, ...
    0.18132296, ...
    0.17810208, ...
    0.17499335, ...
    0.17199103, ...
    0.16908976, ...
    0.16628454, ...
    0.16357067, ...
    0.16094378, ...
    0.15839977, ...
    0.15593477, ...
    0.15354518, ...
    0.15122758, ...
    0.14897878, ...
    0.14679577, ...
    0.14467570, ...
    0.14261590, ...
    0.14061383, ...
    0.13866711, ...
    0.13677347, ...
    0.13493078, ...
    0.13313701, ...
    0.13139025, ...
    0.12968866, ...
    0.12803052, ...
    0.12641419, ...
    0.12483812, ...
    0.12330081, ...
    0.12180086, ...
    0.12033692, ...
    0.11890771, ...
    0.11751201, ...
    0.11614867, ...
    0.11481656, ...
    0.11351463, ...
    0.11224186, ...
    0.11099729, ...
    0.10977999, ...
    0.10858907, ...
    0.10742369, ...
    0.10628304, ...
    0.10516634];



    if (c < 0.1 || c > 10)
        psic = psi(c);
        digammac = psi(1,c);
    else
        i1 = floor(c*10);
        i2 = i1+1;
        psic = t1(i1)+(c-i1/10)*(t1(i2)-t1(i1))/(i2/10-i1/10);
        digammac = t2(i1)+(c-i1/10)*(t2(i2)-t2(i1))/(i2/10-i1/10);

        % or use MATLABs own interpolation methods
        %psic = interp1(t1,c*10,'linear');
        %digammac = interp1(t2,c*10,'linear');
    end
    
    n = length(x);
    onedivn = 1/n;
    wbar = onedivn * sum(log(x));
    tbar = onedivn * sum(x);
    M = log(tbar) - wbar;
    f = log(c)- psic - M;
    g = 1/c - digammac;
end



